library(Signac)
library(Seurat)
library(tidyverse)

# import peak mtx ----------
atac.norm <- Read10X_h5('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/filtered_peak_bc_matrix.h5')
atac.norm[1:5,1:5]

atac.meta <- read_csv('~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/singlecell.csv') |>
  column_to_rownames('barcode')

chrom_assay <- CreateChromatinAssay(
  counts = atac.norm,
  fragments = '~/append-ssd/alaria2/sun2024psoria/ctrl_fullq/outs/fragments.tsv.gz',
  sep = c(":", "-"),
  min.cells = 10,
  min.features = 200
)

aobj <- CreateSeuratObject(
  counts = chrom_assay,
  assay = "peaks",
  meta.data = atac.meta
)

aobj[['peaks']]

granges(aobj)

# keep only standard chromsomes
peaks.keep <- seqnames(granges(aobj)) %in% standardChromosomes(granges(aobj))
aobj <- aobj[as.vector(peaks.keep), ]

# add genome annotation --------
## extract gene annotations from EnsDb packages
# BiocManager::install('EnsDb.Hsapiens.v86') # this is latest hg38 EnsDb
## create from gtf 
# annotations <- GetGRangesFromEnsDb(ensdb = ensdb.v112)

library(AnnotationHub)
ah <- AnnotationHub()

query(ah, "EnsDb.Hsapiens")

ensdb.hs109 <- EnsDb('Homo_sapiens.GRCh38.109.sqlite')

annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86::EnsDb.Hsapiens.v86)

seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg38"

# add the gene information to the object
Annotation(aobj) <- annotations

# QC ------
# compute nucleosome signal score per cell
aobj <- NucleosomeSignal(object = aobj)

# compute TSS enrichment score per cell
aobj <- TSSEnrichment(object = aobj, fast = F)

head(aobj@meta.data)

# compute blacklist frac via Signac 
aobj$blacklist_fraction <- FractionCountsInRegion(
  object = aobj, 
  assay = 'peaks',
  regions = blacklist_hg38_unified
)

# add blacklist ratio and fraction of reads in peaks
aobj$pct_reads_in_peaks <- aobj$peak_region_fragments / aobj$passed_filters * 100
aobj$blacklist_ratio <- aobj$blacklist_region_fragments / aobj$peak_region_fragments

DensityScatter(aobj, x = 'nCount_peaks', y = 'TSS.enrichment', log_x = TRUE, quantiles = TRUE)

aobj$high.tss <- ifelse(aobj$TSS.enrichment > 3, 'High', 'Low')
TSSPlot(aobj, group.by = 'high.tss') + NoLegend()

aobj$nucleosome_group <- ifelse(aobj$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
FragmentHistogram(object = aobj, group.by = 'nucleosome_group') +
  ggtitle('Nucleosome signal QC')

# qc violin plot
VlnPlot(
  object = aobj,
  features = c('nCount_peaks', 'TSS.enrichment', 'blacklist_fraction',
               'nucleosome_signal', 'pct_reads_in_peaks'),
  pt.size = 0,
  ncol = 5
)

# remove cells of bad quality
aobj <- subset(
  x = aobj,
  subset = nCount_peaks > 3000 &
    nCount_peaks < 30000 &
    pct_reads_in_peaks > 15 &
    blacklist_ratio < 0.05 &
    nucleosome_signal < 4 &
    TSS.enrichment > 3
)
aobj

# Normalization & dim reduc & clustering -------
aobj <- RunTFIDF(aobj)
aobj <- FindTopFeatures(aobj, min.cutoff = 'q0')
aobj <- RunSVD(aobj)

# LSI (TFIDF+SVD) component 1 will strongly correlate with seq depth
# So we will exclude it from downstream analysis
DepthCor(aobj)

aobj <- RunUMAP(object = aobj, reduction = 'lsi', dims = 2:30)
aobj <- FindNeighbors(object = aobj, reduction = 'lsi', dims = 2:30)
# SLM clustering
aobj <- FindClusters(object = aobj, verbose = FALSE, algorithm = 3)
aobj |>
  DimPlot(label = TRUE, label.box = T, cols = DiscretePalette(36))

# gene activity matrix ---------
# cost ~5 min
gene.activities <- GeneActivity(aobj)

# add the gene activity matrix to the Seurat object as a new assay and normalize it
aobj[['RNA']] <- CreateAssayObject(counts = gene.activities)
aobj <- NormalizeData(
  object = aobj,
  assay = 'RNA',
  normalization.method = 'LogNormalize',
  scale.factor = median(aobj$nCount_RNA)
)

DefaultAssay(aobj) <- 'RNA'

FeaturePlot(
  object = aobj,
  features = c('KRT1','TRPV3','KRT14','CD3D'),
  pt.size = 0.1,
  max.cutoff = 'q95'
)

skin.marker <- list(
  kera = c('KRT1','KRT10','KRTDAP','KRT14','KRT5','COL17A1'),
  fibroblast = c('FBN1','LUM','COL3A1'),
  macrophage = c('CD163','MRC1','MARCO'),
  T.cell = c('CD3E','CD8A','IL17A','FOXP3'),
  DC = c('CD1A','CD207')
)

aobj |>
  DotPlot(c(skin.marker,'TRPV3'), cluster.idents = T, assay = 'RNA') +
  RotatedAxis()

aobj |>
  DotPlot(c(late.kc,key_cytokine,kegg_mva,'TRPV3'),
          cluster.idents = T,
          assay = 'RNA') +
  RotatedAxis()

aobj |>
  FindAllMarkers(assay = 'RNA',features = c(list_c(skin.marker),'TRPV3'),
                 only.pos = T) |>
  filter(p_val_adj < .05)

library(tidyseurat)

aobj <- aobj |>
  mutate(manual.main = case_when(seurat_clusters == 9 ~ 'T cell',
                                 seurat_clusters %in% c(1,11) ~ 'Macrophage',
                                 seurat_clusters == 10 ~ 'V3h KC',
                                 .default = 'V3l KC'))

# Find differentially accessible peaks ------
# change back to working with peaks instead of gene activities
DefaultAssay(aobj) <- 'peaks'

## simple, acceptable & quick
da_peaks.wilcox <- aobj |>
  FindAllMarkers(only.pos = T, latent.vars = 'nCount_peaks')

head(da_peaks.wilcox)

## use this param will cost >2h for one cluster
da_peaks <- aobj |>
  FindAllMarkers(test.use = 'LR', latent.vars = 'nCount_peaks',
                 only.pos = T)

head(da_peaks)

## identify closest gene from granges
marker.peak.cluster <- da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(p_val_adj < .05) |>
  pull(gene)

marker.peak.close0 <- aobj |> ClosestFeature(marker.peak.cluster)

region.v3 <- marker.peak.close0 |>
  as_tibble() |>
  dplyr::filter(gene_name == 'TRPV3') |>
  pull(query_region)

da_peaks.wilcox |>
  as_tibble() |>
  dplyr::filter(gene == region.v3)

der.v3hvl <- aobj |>
  FindMarkers(ident.1 = 10, ident.2 = 2)

deg.v3hvl <- der.v3hvl |>
  as_tibble(rownames = 'gene') |>
  dplyr::filter(p_val_adj < .05) |>
  pull(gene) |>
  ClosestFeature(object = aobj, regions = _)

deg.v3hvl |>
  as_tibble() |>
  filter(gene_biotype == 'protein_coding') |>
  filter(gene_name %in% c(kegg_mva, late.kc, 'TRPV3'))

# plot genomic region -------
pbmc |> CoveragePlot(
  region = region.v3,
  extend.upstream = 0,
  extend.downstream = 80000,
  assay = 'peaks'
)

hm.skin |> write_rds('mission/FPP/psoriasis/sun24.ctrl.signac.rds')

aobj <- read_rds('mission/FPP/psoriasis/sun24.ctrl.signac.rds')

# transfer anchor from scRNA-seq ----------
sobj_epi <- read_rds('mission/FPP/xiangya_sle_scRNA/data/human_SLE_epi.rds')

transfer.anchors <- FindTransferAnchors(reference = sobj_epi, query = aobj,
                                        reduction = 'cca',query.assay = 'RNA')

predicted.labels <- transfer.anchors |>
  TransferData(refdata = sobj_epi$bill_fine,
               weight.reduction = aobj[['lsi']],
               dims = 2:30)

predicted.labels|> head()

aobj.trans <- AddMetaData(object = aobj, metadata = predicted.labels)

aobj.trans |> DimPlot(group.by = 'predicted.id')

aobj.trans |> DotPlot(c('prediction.score.TRPV3.lo.KC','prediction.score.TRPV3.hi.KC'),
                      scale = F)

aobj.trans |> VlnPlot('prediction.score.max', group.by = 'predicted.id')

# motif analysis ----------
#BiocManager::install(c('motifmatchr','JASPAR2024',
#'TFBSTools','BSgenome.Hsapiens.UCSC.hg38','chromVAR'))
#library(JASPAR2024)
library(TFBSTools)

aobj@meta.data |>
  as_tibble(rownames = '.cell')

pfm <- getMatrixSet(
  x = '~/append-ssd/work/JASPAR2024.sqlite',
  opts = list(collection = "CORE", tax_group = 'vertebrates', all_versions = FALSE)
)

hm.skin <- aobj |>
  AddMotifs(pfm = pfm,
            genome = BSgenome.Hsapiens.UCSC.hg38)

# cost 1-2 min
da_peaks <- hm.skin |> FindMarkers(
  ident.1 = 10,
  ident.2 = 2,
  only.pos = TRUE,
  test.use = 'LR',
  min.pct = 0.05,
  latent.vars = 'nCount_peaks',
  assay = 'peaks'
)

# get top differentially accessible peaks
top.da.peak <- da_peaks |>
  as_tibble(rownames = 'peak') |>
  filter(p_val < .005, pct.1 > .2) |>
  pull(peak)

# find peaks open in KC
open.peaks <- hm.skin |>
  AccessiblePeaks(idents = c(2,10))

# match the overall GC content in the peak set
meta.feature <- hm.skin |>
  GetAssayData(assay = "peaks", layer = "meta.features")
peaks.matched <- MatchRegionStats(
  meta.feature = meta.feature[open.peaks, ],
  query.feature = meta.feature[top.da.peak, ],
  n = 50000
)

# test enrichment
enriched.motifs <- hm.skin |>
  FindMotifs(
  features = top.da.peak,
  assay = 'peaks'
)

enriched.srebp <- enriched.motifs |>
  filter(str_detect(motif.name, 'SREB'), p.adjust < .05) |>
  rownames()

hm.skin |> MotifPlot(motifs = enriched.srebp[1:2])

## Computing motif activities ---------
hm.skin <- hm.skin |> RunChromVAR(genome = BSgenome.Hsapiens.UCSC.hg38)

DefaultAssay(hm.skin) <- 'chromvar'

# look at the activity of SREBF2
hm.skin |> FeaturePlot(
  features = enriched.srebp,
  min.cutoff = 'q10',
  max.cutoff = 'q90'
)

# footprint -------------
## preprocess example data -------------
filepath <- "GSE129785_scATAC-Hematopoiesis-CD34"

peaks <- read_tsv(paste0(filepath, ".peaks.txt.gz"))
cells <- read_tsv(paste0(filepath, ".cell_barcodes.txt.gz"))
cells <- cells |>
  column_to_rownames('Barcodes')

mtx <- Matrix::readMM(file = paste0(filepath, ".mtx.gz"))
mtx <- as(object = mtx, Class = "CsparseMatrix")
colnames(mtx) <- rownames(cells)
rownames(mtx) <- peaks$Feature

bone_assay <- CreateChromatinAssay(
  counts = mtx,
  min.cells = 5,
  fragments = "GSM3722029_CD34_Progenitors_Rep1_fragments.tsv.gz",
  sep = c("_", "_")
)

mtx[1:5,1:5]

cells |> rownames() |> head()

bone_assay@counts |>
  glimpse()

bone <- CreateSeuratObject(
  counts = bone_assay,
  meta.data = cells,
  assay = "ATAC"
)

bone@meta.data |>
  head()

bone$Group |> table()

# The dataset contains multiple cell types
# We can subset to include just one replicate of CD34+ progenitor cells
bone <- bone[, bone$Group == "CD34_Progenitors_Rep1"]

# add cell type annotations from the original paper
cluster_names <- c("HSC",   "MEP",  "CMP-BMP",  "LMPP", "CLP",  "Pro-B",    "Pre-B",    "GMP",
                   "MDP",    "pDC",  "cDC",  "Monocyte-1",   "Monocyte-2",   "Naive-B",  "Memory-B",
                   "Plasma-cell",    "Basophil", "Immature-NK",  "Mature-NK1",   "Mature-NK2",   "Naive-CD4-T1",
                   "Naive-CD4-T2",   "Naive-Treg",   "Memory-CD4-T", "Treg", "Naive-CD8-T1", "Naive-CD8-T2",
                   "Naive-CD8-T3",   "Central-memory-CD8-T", "Effector-memory-CD8-T",    "Gamma delta T")
num.labels <- length(cluster_names)
names(cluster_names) <- paste0( rep("Cluster", num.labels), seq(num.labels) )
new.md <- cluster_names[as.character(bone$Clusters)]
names(new.md) <- Cells(bone)
bone$celltype <- new.md

bone[["ATAC"]]

# latest hg19 ensdb
library(EnsDb.Hsapiens.v75)

# extract gene annotations from EnsDb
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)

# change to UCSC style since the data was mapped to hg19
seqlevels(annotations) <- paste0('chr', seqlevels(annotations))
genome(annotations) <- "hg19"

# add the gene information to the object
Annotation(bone) <- annotations

bone <- TSSEnrichment(bone)
bone <- NucleosomeSignal(bone)
bone$blacklist_fraction <- FractionCountsInRegion(bone, regions = blacklist_hg19)

bone <- bone[, (bone$nCount_ATAC < 50000) &
               (bone$TSS.enrichment > 2) & 
               (bone$nucleosome_signal < 5)]

bone <- FindTopFeatures(bone, min.cells = 10)
bone <- RunTFIDF(bone)
bone <- RunSVD(bone, n = 100)
DepthCor(bone)
bone <- RunUMAP(
  bone,
  reduction = "lsi",
  dims = 2:50,
  reduction.name = "UMAP"
)
bone <- FindNeighbors(bone, dims = 2:50, reduction = "lsi")
bone <- FindClusters(bone, resolution = 0.8, algorithm = 3)
DimPlot(bone, label = TRUE) + NoLegend()

for(i in levels(bone)) {
  cells_to_reid <- WhichCells(bone, idents = i)
  newid <- names(sort(table(bone$celltype[cells_to_reid]),decreasing=TRUE))[1]
  Idents(bone, cells = cells_to_reid) <- newid
}

DimPlot(bone, label = TRUE)

DefaultAssay(bone) <- "ATAC"

erythroid <- bone[,  bone$assigned_celltype %in% c("HSC", "MEP", "CMP-BMP")]
lymphoid <- bone[, bone$assigned_celltype %in% c("HSC", "LMPP", "GMP", "CLP", "Pro-B", "pDC", "MDP", "GMP")]

## footprinting ---------
library(motifmatchr)
library(JASPAR2020)
library(TFBSTools)
library(BSgenome.Hsapiens.UCSC.hg19)

# extract position frequency matrices for the motifs
pwm <- getMatrixSet(
  x = JASPAR2020,
  opts = list(species = 9606, all_versions = FALSE)
)

# add motif information
bone <- AddMotifs(bone, genome = BSgenome.Hsapiens.UCSC.hg19, pfm = pfm)

# gather the footprinting information for sets of motifs
bone$celltype <- Idents(bone)

bone.mep.hsc <- bone |>
  dplyr::filter(celltype %in% c('MEP','HSC'))

Idents(bone.mep.hsc) |> table()

bone.mep.hsc <- Footprint(
  object = bone.mep.hsc,
  motif.name = c("GATA2"),
  in.peaks = T, upstream = 80, downstream = 80,
  genome = BSgenome.Hsapiens.UCSC.hg19
)

# plot the footprint data for each group of cells
p2 <- PlotFootprint(bone.mep.hsc, features = c("GATA2"))

p2 + patchwork::plot_layout(ncol = 1)

p2[[1]]

bone$celltype |> table()

## check threshold for deviation? ----------
bone <- bone |>
  RunChromVAR(genome = BSgenome.Hsapiens.UCSC.hg19)

chrvr.mep.hsc <- bone |>
  FindMarkers(
    ident.1 = 'MEP',
    ident.2 = 'CMP-BMP',
    only.pos = F,
    mean.fxn = rowMeans,
    fc.name = "avg_diff",
    assay = 'chromvar'
  )

motif.meta <- bone |>
  GetMotifData(slot = 'motif.names') |>
  as_tibble() |>
  pivot_longer(everything())

motif.gata2 <- motif.meta |>
  dplyr::filter(str_detect(value, 'GATA2'))

motif.gata2

chrvr.mep.hsc |>
  as_tibble(rownames = 'id') |>
  dplyr::filter(id == motif.gata2$name)
